Mon. Not. R. Astron. Soc. 000, ITl- Il3l (00001 Printed 3 March 2013 (MN WT^i style file v2.2) 



Pseudo-bulge formation via major mergers 



(N 

o 



< 

o 

43 
6 

o3 



> 
o 

(N 



13 



J. A. Keselman 1 *, and A. Nusser 1,2 

1 Physics department, Technion, Haifa 32000, Israel 

2 Asher Space Research Institute, Technion, Haifa 32000, Israel 



3 March 2013 



ABSTRACT 

It is widely accepted that within the framework of ACDM a significant fraction of 
giant-disk galaxies has recently experienced a violent galactic merger. We present 
numerical simulations of such major mergers of gas-rich pure disk galaxies, and focus 
on the innermost stellar component (bulge) of the disk remnants. The simulations 
have high spatial and mass resolutions, and resolve regions deep enough to allow 
bulge classification according to standard kinematical and structural characteristics. 
In agreement with recent studies we find that these bulges are dominated by stars 
formed in the final coalescence process. In contrast to the common interpretation of 
such components as classical bulges (i.e. similar to intermediate luminosity ellipticals), 
we find they are supported by highly coherent rotations and have Sersic indices n < 2, 
a result leading to their classification as pseudo-bulges. Pseudo-bulge formation by 
gas rich major mergers of pure disks is a novel mode of pseudo-bulge formation; It 
complements pseudo-bulge growth by secular evolution, and it could help explain the 
high fractions of classically bulge-less giant disk galaxies, and pseudo-bulges found in 
giant Sc galaxies. 
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1 INTRODUCTION 

A galactic bulge in a spiral galaxy is morphologically de- 
fined as the excess of light in the central region of the galac- 
tic disk. With this simple definition, recent numerical sim- 
ulations which include sub-grid Inter-Stellar Medium (ISM) 
physics are beginni ng to successfully reproduce the observed 
bulge masses (e.g. [S teirmictz' fc Navarro] 120021 ; iParrv et alj 
120091 ; [Hopkins et al.ll2O10f ). However, bulges have complex 
dynamical properties that are still inaccessible by means of 
current simulations. Detailed observations of nearby galax- 
ies reveal spiral galaxies harbouring classical bulges which 
are supported by random motions of their stellar compo- 
nents. Others spirals, however, show evidence for pseudo- 
bulges, those are bulges that are mainly supported by co- 
herent rotation. Current models for bulge formation remain 
challenged when confronted with the observed division into 
classical and pseudo-bulges. The key issue is that there are 
fewer classical bulges than expected from models for bulge 
formation by major mergers in the co ntext of the standard 
hierarchical ACDM paradigm (see e.g. lWhite fe Reeslll978l ; 
ICole et al ]|l994|Peebleslll993h . In these models major merg- 
ing are ubiquitous and are thought to prod uce remnants that 
are ce ntrally dominated by classical bulges. [Kormendv et alj 
studied galaxies with Vdrc > 150 km s 1 within a 
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sphere of radius 8Mpc centered on our Galaxy. They find 
that 11 out of 1 9 gala xies show no sign of classical bulges. 
iFisher fc Drorvl (|201ll ). using Spitzer 3.6pm detector and 
Hubble space telescope data, find that within an UMpc 
sphere the dominant galaxy type has pure disk character- 
istics, whether counting by number, star formation rate, or 
stellar mass. 

The classification into classical and pseudo-bulges could 
be ma de on the basis of a variety of observational charac ter- 
istics jKormendv fc Kennicuttll2004lAthanassoulall2005t ). In 
addition to kinematics, of particular importance is the distri- 
bution of stars in the galactic components. The characterisa- 
tio n of the stellar dis tribution is assesed by the Sersic index 
n (|Sersidll963l .l 1968), used in the following parametrization 



S(r) = I e exp 



1/r, 



- 1 



(1) 



of the the observed (surface) stellar intensity profile, E(r). 
This form includes two free parameters, R e is an effective 
radius enclosing half the modelled mass, and J e which is the 
intensity at R e - The parameter b n is defined by the latter 
requirement, i.e. it is calculated from I e and R e . The stellar 
distribution in both bulges and pseudo-bulges is well ap- 
proximated by the form[TJ typically with n > 3 for classical 
bulges, and n < 2 for pseudo-bulges, closer to exponential 
disks (n = 1). 

Here we explore the possibility that violent major merg- 
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ing could actually lead to the formation of pseudo-bulges. 
We find that gas rich major mergers of pure disk galaxies 
could indeed produce inner cores with morphological and 
kinematical properties of the observed pseudo-bulges. This 
does not mean that violent mergers will be the main mode 
of pseudo-bulge formation since it requires the existence of 
a significant population of pure disk galaxies. However, it 
could greatly mitigate the tension between standard bulge 
formation models and observations. The physical conditions 
for pseudo-bulge formation by secular growth during the qui- 
escent later phase of galaxy evolution, remain more likely to 
be satisfied. 

Like most recent studies of galactic bulge formation 
via mergers, we take into account key physical processes 
such as radiative dissipation, star formation, and feedback. 
The reason these are needed is that classical bulges ob- 
tained from dry major mergers do not resemble interme- 
diate luminosity ellipticals like the observed bulges do, 
i.e. flattening by rotation (jDavis et al.l Il983|). Sersic pa- 
rame ter n that scales with lumin osity (I Graham fc Guzman! 
120031 ). color-magnitude relations (jBalcells fc Peletierill994 
and metallicity-luminosity relations (jjablonka et al.l 1 1996 
Bulges created by dry mergers are generally more similar to 
giant ellipticals or diffuse s tellar halos. This can be seen, for 
e.g., in Fig. 5 in ICox et al.l(|2006l l. where remnants of gas-less 
galaxies show no sign of a significant inner component. 

We now review a few recent studies of galactic bulge for- 
mation via mergers, and focus on difficulties reg arding the 
classifi cation into classical versus pseudo-bulges. ICox et all 
l|2006f ) study remnants of galactic mergers and analyse them 
using observational methods, while treating the whole rem- 
nant as a single component. The whole remnant being a 
single component, they find large effective radii with mean 
of 4kpc, a number significantly larger than measured in 
this work. Also, their simulations do not have the required 
spatial resol ution to measure the Sersic index of the inner 
co mponent. Hopkins et al.l (|2010h uses simulations similar 
to lCox et al. (|2006l ) . In their decomposition of the remnants 
into various components, bulges are assigned only those par- 
ticles with weak rotational support. Essentially, any disk-like 
bulge is not trea t ed, an d is assumed to be part of the disk. 
iRobertson et al.l l)2006h use similar simulations, but adopt 
two methods to decompose the remnants. First, much like 
what we do here, a remnant is decomposed into a bulge, 
thin disk, thick disk, and a spheroid, according to their stel- 
lar surface density. This decomposition is then used mainly 
to derive relative masses. Second, they decompose the rem- 
nant into old stars, that formed a few hundred Myrs before 
the peak of the star formation rate (hereafter SFR), and 
young stars, formed a few hundred Myrs afterwards. They 
find that young stars are highly rotationally supported (disk 
like) while old stars are not (classical bulge). The composi- 
tion of the remnant core (whether it is composed mainly 
of young or old stars) is not addressed. The stellar popula- 
tion formed during the main SFR peak is not analysed as a 
separate component. However, as we show here, these stars 
may have a maj or contribution t o the central density of the 
galaxy (see also lCox et al.ll2006t ). constituting a significant 
component which could be classified as a pseudo-, rather 
than a classical bulge. Pure bulge-dominated spheroids, with 
rotating central components formed via dissipative inflows 
and "in-situ" star formation by the inflowing gas in the fi- 
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None of these studies directly address the question of 
whether pseudo-bulges could form via major mergers. Here, 
we revisit the problem of bulge formation by major mergers 
through simulations, including the usual essential processes 
of star formation, SN feedback and radiative dissipation of 
gas. Our simulations are of sufficiently high spatial and mass 
resolution to allow a realistic analysis of the simulated stellar 
cores in more details than done before. We find that despite 
their violence, dissipative major mergers can produce a cen- 
tral stellar component with highly coherent rotation. 

The structure of the paper is as follows. In ij2]we give 
a high level description of the numerical experiment. In Ej3] 
we describe the numerical model: we start with a brief de- 
scription of the simulation code, then in i)3.1l we describe 
the galactic model used to create the simulated galaxies. 
In ij4] we describe the simulations (orbits, initial positions, 
integration time, and so on), in ij5]the structural and kine- 
matic properties of the simulated remnants, in ij6] numerical 
issues, and we conclude with a general discussion about the 
implications of our findings in jJT] We generate galaxies with 
massive stellar and gaseous disks in equilibrium using a new 
method, described in detail in appendix 1X1 



2 THE GENERAL SET-UP 

We wish to derive conditions under which merging of pure 
disk galaxies is likely to form pseudo-bulges at the cores 
of disk-like galaxies. We address this problem using a suite 
of merging simulations of initially isolated multi-component 
galaxies which include gas, stars and collisionless dark 
haloes. The wide range of possible configurations of the 
merging galaxies is ideally probed in a large scale cosmolog- 
ical simulation. Unfortunately such a simulation including 
all relevant physical effects, and the sub-kpc resolution re- 
quired to model the cores of galaxies, is unavailable. Instead 
we aim at extracting robust general conclusions from a few 
detailed merging simulations. 

We perform seven simulations of mergers of identical 
pairs of pure disk galaxies. Prior to the merging, the orig- 
inal galaxies fall into four configurations corresponding to 
two choices for the gas mass fraction in the disk and two 
choices for the effective equation of state of the ISM. These 
original galaxies have been generated in isolation, all shar- 
ing the same choices for the masses in the baryonic and 
dark components, the halo density profile, and the disk scale 
lengths (vertical and radial) . They also have the same num- 
ber of baryonic and DM simulation particles, and the same 
criteria for setting the force softening lengths. The original 
galaxies all have r2oo = 232 kpc, roughly like the Milky- Way 
ijYin et al.ll2009h . As mentioned in bli the DM distribution 
does not strictly follow an NFW profile, however, at small 
radii, the profile is close to an NFW form with c NFW = 9. 
The values of all relevant parameters are listed in table Q] 
and explained in details in the caption. The table also lists 
a fifth galaxy (gO) which is only used to test the validity of 
the scheme used in generating equilibrium configurations of 
collisionless stellar disk embedded in a spherical DM halo. 
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Figure 1. The star formation rate as a function of time since the beginning of the simulation. 



3 NUMERICAL MODEL 



3.1 Galactic model 



We use a modifie d version of t he Gadget2 hydrodynami- 
cal N-body code (|Springell [20051 ). which includes radiative 
dissipation of gas energy, star formation an d SN feedback, 



Springel fc HernauistJ l|2003h . In this model, 



as described in 

The ISM is cooled by line emissions of a primordial gaseous 
mix, and heated by SN energetic feedback. It can be shown 
that under some reasonable assumptions, in star-forming re- 
gions, the ISM is in a stable phase, and exhibits an effective 
equation of state (EEOS), in which the pressure is a func- 
tion of density. This is the state of the ISM assumed in the 
simulations, in regions with gas density above a threshold 
pth - Since this is a sub-grid model, and the feedback mecha- 
nism is not completely understood, there is some freedom in 
choosing the ISM pressurization. The parameter Q is used 
to linearly interpolate the EEOS between the full model and 
an isothermal EOS (Q=0 means there is no feedback, and 
Q — 1 means the fu ll model is used, wit h the chosen ef- 
ficiency parameters). iHopkins fc Quataertl (|201Cf ) compare 
various choices for the Q parameter to observations, with the 
conclusion that Q ~ 0.1-0.3 pr ovide the best match. We use 
the same parameters as used bv lSpringel et al l (|2005h . where 
the efficiency of the evaporation process Aq = 4000, the mass 
fraction of stars that are short-lived and die as supernovae 
P — 0.1, the cold gas cloud temperature T c i ou d — 1000K, the 
effective SN temperature Tsn = 4 x 10 8 K, and the gas con- 
sumption time-scale t* = 8.4 Gyrs. This choice gives a stellar 
formation rate of ~ 1 Mq per year for isolated galaxies with 
structure similar to the Milky- Way. In the numerical repre- 
sentation, each gaseous particle may spawn up to n 3 stellar 
particles of mass m g /n 3 where m g is the initial mass of the 
gas particle. In order to assess the numerical effect of n s We 
use simulations both with n s — 6, and n s — 2. Throughout 
this study we adopt the Hubble parameter h — 0.7. 



We follow Springel ct al. (20 051) and model the dark matter 
(DM) halo as a lHernauisti(|l99Gi ) profile 



p(r) = 



M d 



2n r(r + a)' 



(2) 



where Mdm is the halo mass for r — > oo, and a is a scale 
length. The length scale a can be matched to the concentra- 
tion parameter c and the scale r„ in the Navarro, Frenk fc 
White density profile, p NFW (r) |Navarro et al.lll997l , Il996bl . 
here after NFW) by demanding that p(r) = p NPW (r) for 
r <C J"2oo, where r2oo is the radius within which the mean 
DM density is 200 times the critical density of the universe. 
This yields the relation 



nV2[m(l + c) - c/(l + c) 



(3) 



Here we use for all galaxies a = 42kpc, r s — 25kpc, and 
c=9. 

We adopt an exponential stellar disk with a length scale 
h, and with vertical structure of an isothermal sheet with a 
constant height scale zo 



P*(R, z) 



4,-kzK 2 



sech 



2 (—) 

\2z a J 



exp 



(4) 



where R is the position vector projected in the disk plane. 
The gaseous component has the same radial structure as the 
stellar component, and its vertical structure is determined 
by the self-consistent solution of the hydrostatic equilibrium 
and Poison equations; We generate the galaxies with both 
collisional and collisionless components in equilibrium using 
a new method, as explained in appendix [X] The used values 
of all relevant parameters are listed in table [1] and explained 
in details in the caption. Note that some parameter choices, 
like the disk gas fractions, are similar to those used in 
previ ous work (|Springel fc Hernqui st 2005; Robertson et alj 
2006), and are essential to achieve disk-like remnants. 



J. A. Keselman & A. Nusser 




g70q1!o3n2 
= 2.0£yr 




DM | g70q1i)1 
t= 2.0<Sy 




g70q1|o2 
t=2.0<3y 




-1.5 -1 -0.5 0.5 1 1.5 2 
Log(r) [kpc/h] 



-1.5 -1 -0.5 0.5 1 1.5 2 -2 -1.5 -1 -0.5 0.5 1 1.5 2 -2 -1.5 -1 -0.5 0.5 1 1.5 2 
Log(r) [kpc/h] Log(r) [kpc/h] Log(r) [kpc/h] 



350 

300 

250 

3» 200 
E 

7 150 

100 

50 

350 

300 

250 

J» 200 
E 

~ 150 
100 
50 




g3|0q1o3 
t=|1.7Gyr 


g3pq025o3 
t={1.7Gyr 


g7|0q1o3 
t=|2.0Gyr 


g70q025o3 
t = ;1.7Gyr 


























^ \ — h- 1 1 1 












g70q1o3n2 
t = £.0Gyr 


g70q1o1 
nfct^ p.OGyr 


i g7pq1o2 
U=£.0Gyr 


0.2 0.4 0.6 0.8 1.0 1.2 
r [kpc/h] 




- J 




I 


^ — "N. / 






Ijf 











0.2 0.4 0.6 0.8 1 1.2 0.2 0.4 0.6 0.8 1 1.2 0.2 0.4 0.6 0.8 1 1.2 
r [kpc/h] r [kpc/h] r [kpc/h] 



Figure 2. Equilibrium of the remnants in terms of mass profile and velocity moments. Top panel: Cumulative mass within radius r 



from the center. Solid lines are measured at t 



and dashed lines are measured 143 million years later. Here, G is for gas, SD for 



stellar disk, SF for stars formed during the collision, and DM for dark-matter. Vertical dashed lines denote formal numerical resolution. 
Bottom panel: Velocity moments of remnants. Velocities are measured like an observer along the edge-on line of sight. Black lines denote 
rotational velocity, and red lines central velocity dispersions. Vertical dashed lines denote formal numerical resolution. 



4 THE MERGING SIMULATIONS 

We perform seven merging simulations, each of identical 
pairs of pure disk galaxies, as described in the previous sec- 
tion. We use three distinct orbits. Initially, galaxies are sep- 
arated by 143 kpc, each galactic disk in the merging pair is 
oriented in the same direction, defining the y axis, and initial 
velocities are as summarised in table [21 Initially all galaxies 
rotate in the same direction, and their angular momentum 
is negatively oriented in the y axis. 



The galaxies begin to rotate around the center of mass, 
while tidal forces and dynamical friction bring them closer 
with time. When galaxies are close enough, their galactic 
disks collide for the first time, and a first starburst is trig- 
gered. Close encounters of the disks significantly enhance the 
gas density, setting ripe conditions for starbursts to occur. 
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Table 1. Parameters used in the generation of the original galaxies. N^ m : number of DM particles. TV*: number of stellar particles in 
the disk. N g : number of gaseous particles in the disk. N ap h. number of neighbours used in the SPH calculation, f g : gas mass fraction in 
the disk, q: pressurization parameter (has some effect on the effective equation of state, see S|3). a.: scale length of the DM halo, h: radial 
scale length of the disk, zo: height of the disk. Mdm'- DM mass (in 10 10 Mq). Ma: total (gas+stars) disk mass (in 10 10 Mq). £dm- force 
softening length for the DM particles, e*: force softening for the initial stellar particles. e s f. force softening for the newly formed stellar 
particles. e g : gravitational force softening for the gas particles, equal to min(e ap h, 0.07 kpc) where e ap h is determined according to the 
nearest 64 particles. 
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Table 2. Parameters used in the merging simulations. x,y & z: initial relative positions between the two merging galaxies. v x ,v y & v z : 
initial relative velocities. n s : the maximum number of stellar particles that a gas particle may spawn. t r f. is an approximate estimate 
for the formation time of the remnant. t max : the time at which simulation was stopped. 



This behaviour is shown in Fig. [TJ where most simulations 
show a secondary starburst peak before a main one. Most 
of the stellar mass of the core of the remnant is produced in 
the main starburst corresponding to the final coalescence of 
the disks. 

We denote the approximate time when the remnant is 
formed by t r f, and loosely define it as the time when the 
core mass and velocity structure is no longer changing sig- 
nificantly during a dynamical time (see Fig. [2} We continue 
the simulations up to time tj. These times are given in ta- 
ble Fig. [3] showcases a remnant at t/. This remnant con- 
tains a thick disk made exclusively of newly formed stars. 
This behaviour has been see n in numerous othe r stud ies, e.g. 
ISpringel fc HernquistJ (|2005T ); iRobertson et al.l (|200rj ). 

In the analysis of the remnants, given below, we only 
consider separations larger than 4e, where e is the spline 
smoothing of the gravitational force. For this thresholdthe 
simulated remnants have a relaxation time longer than the 
simulated time, see Fig. [3] and the results are stable with 
time. Further, at this separation the force between two sim- 
ulated particles is exactly Newtonian (in Gadget this hap- 
pens at 2.5e). Repeating the analysis with a threshold at 
3e, yields similar results. See Sj6]for a detailed discussion of 
numerical considerations. 



5 STRUCTURE OF THE REMNANTS 

The remnants contain a compact pseudo-bulge, a thick stel- 
lar disk, a large gaseous disk, and an extended diffuse stellar 
halo consisting of old stellar material originating from pro- 
genitor disk s. This structure is typical o f such simulations, 
see fo r e.g. ISpringel fc Hernquistl (|2005l ); iRobertson et al.l 
(|200rJ ). A thick stellar disk is showcased in Fig. [3] ac- 
cording to Fig. [5] (top panel) the stellar halo consist of 
~ 1.5 x 1O 1O M in regions outside the stellar disk 6 < r < 
60 kpc, and compared to the disk like components it is dif- 
fuse, and it is similar to known diffuse stellar haloes of disk 
galaxies in both mass and extent- the Milky Way stellar halo 
mass b eing ~ 1 x 1 9 and extended up to ~ 40 kpc see 
for e.g. iTumlinsonl (|2010l ). 

Our analysis of the simulated remnants aims at mim- 
icking the actual observations in as many details as possi- 
ble. Thus, for each remnant, two-dimensional (2D) images 
are obtained projecting the three-dimensional (3D) parti- 
cle positions and velocities onto 1000 planes corresponding, 
respectively, to random viewing directions limited by a max- 
imum inclination of n/6 as to preserve a significant edge-on 
velocity component. On a random galaxy sample this means 
we only measure 1 out of 3 galaxies. 
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Figure 3. Top panel: Contour plot of log the stellar surface 
density (in units of 10 10 M Q /i/kpc 2 ) for g70qlo3 at t = 4.0Gyrs. 
Bottom panel: Line of sight velocity contours of g70qlo3 at t = 
4.0Gyrs. In both panels the galaxy is seen edge-on. 



In order to determine the structural parameters R e and 
n, as appearing in Eq.[T] we decompose the projected density 
of all stellar matter in the remnants as a sum of two circu- 
lar Sersic functions. The fit is performed up to the formal 
simulation spatial resolution. A decomposition of remnants 
into two Sersic functions is shown in Fig. 

We measure the ellipticity of the core. We do so by 
matching an ellipse to the isodensity contour containing 
most of the core mass. We define this mass as the projected 
mass within the circle with radius Rmax, which in turn is 
defined as the circular radius where the core density is 5% 
stronger than the outer component. This is a somewhat arbi- 
trary choice, but we find results to be only weakly dependent 
on it. The results are shown in Fig. [S] 

Observers usually fit images with a 2D decomposition, 




- - - 0.0 0.5 1.0 1.5 

Log(r) [kpc/h] 

Figure 4. Top panel: Relaxation time of newly formed stellar 
components in at tt, as function of distance from the center. 
Dashed black curve is for g30qfo3, dotted black for g70qlo3n2, 
blue for g30q025o3, red for g70q025o3, cyan for g70qfol, pur- 
ple for g70qfo2, and green for g70qfo3. The vertical solid black 
line denotes r = 4e s j-, which we consider as the spatial resolu- 
tion limit. Bottom panel: Dynamical time of newly formed stellar 
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rather than a circular one. In such decomposition each com 
ponent may be described by the following Sersic function: 

— \ I/" 

Ti 



E(r-) = I e exp 



(-^y* 2 + (a - ^)2/) 2 ) 



i 



(5) 



where the main remnant axis is defined as the x axis. This 
decomposition gives different Sersic parameters when com- 
pared to the circular one. We calculate the Sersic parameters 
of the two dimensional decomposition using the following 
method: we find the structural parameters of the observed 
image, the ellipticity e, and the circular Sersic Rl D , and n 
Next we construct a 2D particle realization using the pre- 
viously measured ellipticity e, and random two-dimensional 
Sersic parameters R 2D , and n 2D . The purpose is to find such 
2D parameters that will reproduce the observed image. Thus 
we build many such 2D particle realizations, until we find 
the one that reproduces the circular R\ D , and n 1D . Two- 
dimensional parameters of the simulated remnants are given 
in table |3] In general, we find that n 2D , for n 2D ~ 1, is usu- 
ally larger than n 1D for the same galaxy. This can be seen 
in Fig. [7] In that figure, the over-prediction of a parameter, 
say n, is defined as 



op% 



x 100% 



(6) 



where n c , and U2d are the circular and the two-dimensional 
Sersic parameters of the remnants. 

The most important structural feature of the remnants 
is the small Sersic index, n < 1, of the core. Note that the 
core created in simulation g30qlo3 is very small, and is not 
dominant over the outer galactic component within the sim- 
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Figure 5. Density profile of the remnants as a function of distance from the center. Green dotted lines are the Sersic fit. Vertical dashed 
lines denote formal numerical resolution, and red cross is Umax- 



simulation 


e 


Re 


n 


C/T 


name 




kpc 






g70qlol 


0.1 


0.4 


0.79 


0.41 


q70qlo3 


0.4 


0.29 


0.57 


0.25 


g70qlo3n2 


0.45 


0.37 


0.33 


0.27 


g70qlo2 


0.4 


0.26 


0.67 


0.42 


g70q025o3 


0.14 


0.11 


0.85 


0.1 


g30q025o3 


0.33 


0.1 


0.79 


0.03 



Table 3. Two-dimensional Sersic parameters of the remnants 
core at tf. The last column give the core to total light ratio. 



significant fraction of viewing angles. This is true at least at 
a certain period of time in the galactic evolution. 

Fig. [8] is plotted using an observationally inspired 
method: Velocity is measured along a thin slit aligned with 
the remnant image main axis. The slit ratio is 12 : 1, and 
length 2R ma x. It is divided into bins of 3000 particles. We 
measure Vmax as the maximal velocity along the slit, and a 
as the velocity dispersion in the inner square of side length 
Rmax/6- It is clear that a significant fraction of the viewing 
angles leads to a classification of the core as pseudo-bulge. 
As further analyses, we plot the fraction of viewing angles 
that leads to such classification as function of time. See Fig. 
[5] Even for our worst-case merger, and after 4Gyrs from the 
start of the simulation, more than 20% of the viewing angles 
lead to a pseudo-bulge classification. 



ulation spatial resolution. Thus, we do not include this core 
in most of the structural analyses. 

We quantify the kinematic structure of the rem- 
nants by positioning the m in the anisotropy diagram 
(|Binnev fc Tremaind Il987l ) . A point in this diagram con- 
sists of the ratio a /Vmax versus the ellipticity e, where Vmax 
is the maximal rotational velocity of the central component. 
The information given by such a point may be interpreted 
as a measure of how much the gravitational system under 
inspection is supported by rotation. The solid lines in Fig. 
[8] under some simplifying assumptions, are expected values 
of a /Vmax for highly rotationally supported systems with 
ellipticity e. Systems above this line are considered to be 
rotationally supported, and systems below it- supported by 
random motions. In Fig. [8] we plot the different simulated 
remnants as viewed from 1000 different angles, limited by 
inclination as described before. This plot shows that almost 
all simulated remnants appear rotationally supported for a 



6 NUMERICAL CONSIDERATIONS 

In Fig. [2] we show that the newly formed cores are indeed 
close to equili brium, with the velocity mome nts in the ob- 
served range (jKormendv fc Illingworthlll982l ). Here we fur- 
ther investigate possible numerical effects related to the fi- 
nite relaxation time and smoothing. The dynamical time 
scales of the remnants at their formation times i r / are shown 
in Fig. |3] The dynamical time at distance r from the center 
is calculated as td(r) = R 1/2 /a 1/2 , where R 1/2 is the radius 
enclosing half the mass within r, and a 1/2 is the velocity 
dispersion within R 1/2 - The remnants core relaxation time 
at time tf and near the center (up to the resolution limit) 
is longer tha n the simulated time. Th e relaxation time is 
calculated as (|Farouki fc Salpeterlll994 ) 

N 
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Figure 6. Sersic parameters of the remnants. Dotted red, 
solid black, and green dashed curves, correspond to g70qlo3, 
g70qlo3n2, and g70qlol, respectivly. 



where N is the number of particles within r, F is an empir- 
ical dimensionless factor, td is the dynamical time at r, and 
A = R 1/2 /e. We use F = 8. 

Two Numerical considerations in the simulations pre- 
sented in this paper may be considered as marginal; First, 
the relaxation time. Fig. [4] and Table [2] show the relaxation 
time of the remnants at the final time of the simulation, 
tf. At the innermost studied region r = 4e, the relaxation 
time, t r , is always larger than the simulated time tf. While 
t r ~ 3t/ for simulation g70q025o3, it is close to the simu- 
lated time for g30qlo3, and it is in between these values for 
the other simulations. Of course tf is longer than the age of 
the remnant which formed about 1 Gyrs after the beginning 
of the simulation. Second, we analyse results in a region 1.6 
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Figure 7. Systematic over-predictions of Sersic parameters as 
a function of the remnants ellipticity, measured as described in 
the text. Dashed red and solid lines correspond to R e , and n 
respectively. Heavy and light lines are for R m i n = 1.1 X R e , and 
Rmin = 0. Top panel is for n = 0.7, and bottom for n = 3.5. 



times larger than the distance at which gravitational force 
in the simulation becomes exactly Newtonian. 

We now discuss in more det ail the above considerations . 
The gravothermal catastrophe (|Lvnden-Bell fc Woodlll96ct ) 
is one of the most important gravitational effects driven by 
collisions, either physical or numerical. Due to the negative 
heat capacity of such systems, if the center is hotter than 
the outskirts, heat flows outwords resulting in further heat- 
ing and contraction of the central region. As the core col- 
lapses, it also r otates faster even though angular momentum 
flows outwards |Lvnden-Bell fc Pri nglc 1974; Hachisu 19791 ). 
The magnitude of these effects after a si ngle relaxati o n tim e 
is modest, as can be seen in Fig. 4 in lErnst et al.l (|2007h . 
ISpurzem et"aH l|2003h studied the evolution of the ratio 
v rot I a and fo und that it monotonic ally decreases with time. 
Although this lSpurzem et all (|200St ) result seems favourable 
to ours, in order to classify a bulge as pseudo w e need 
to know its ellipticity. Fig. 7a in lEinsel fc Spurzeml (|l999T ) 
shows how the dynamical ellipticity = yjl — (1 — e) 2 
evolves with time. The evolution of the cores with maxi- 
mal mean ellipticity (i.e. ellipticity averaged over viewing 
angles), e ~ 0.5, within a single relaxation time, is mod- 
est. To further assess the two-body relaxation effects on our 
results we use compare simulation g70qlo3n2 and g70qlo3. 
These two simulations are identical in every respect, except 
g70qlo3n2 contains a third of the newly formed stellar par- 
ticles. Fig. [9] shows that the fraction of viewing angles lead- 
ing to a pseudo-bulge classification of the core of simulation 
g70qlo3n2 is lower than that in simulation g70qlo3. Thus 
we conclude that two-body relaxation effects, when using 
the anisotropy diagram, are more likely to lead to classi- 
cal bulges rather than pseudo. Hence our identification of 
pseudo-bulges is robust. 
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Figure 8. Anisotropy diagram of remnants core. Each galaxy is viewed from 1000 random angles, each viewing angle is represented by 
a point. 



Two numerical effects, the gravitational force smooth- 
ing, and core collapse due to relaxation, compete to decrease 
and increase, respectively, Sersic n. Under these circum- 
stances, we would like to understand whether or not the 
density profile can be considered as physically converged in 
the region r < 4e. Unfortunately this is a grey area, as for 
our case there is no goo d agreement between different stud- 
ies: iMoore et alj (|l998l ) claim that the minimum resolved 
scales should be no larger than half the mean interparticle 
separation within the virial radius, and that several thou- 
sands of particles are needed to resolve the innermost re- 
gions of collisionless systems. T hese criteria are all met in 
the simulations here. However, iGhigna et ail (|200Ch claim 
that, in addition to a large number of particles, convergence 
requires that, in addition to a large number of particles, the 
minimal resolved radius is at least as large as three times 
the distance at which gravitational forces between two par- 
ticles becomes Newtonian. This means we should have used 
at least r > 7e, but these regions already miss some of the 
most interesting features of the bulge. However, we did not 
see any difference when we anal ysed the results with r > 3e 
compared to r > 4e. The study of lKlypin et al.l (|200ll ) claims 
that regions containing 200 particles can alread y be reliably 
resolved provided r > 4e. In lPower et alj l|2003l ) a new con- 
dition for convergence is presented. This condition is based 
on a maximum characteristic acceleration imprinted by soft- 
ening. It is empirically studied for NFW haloes, and it is not 
clear how it is applicable to a galactic bulge. 

As a final note, the bulges presented in this study are 
too bright when compared to observed bulges. They are out- 
liers (from above) in the Kormend y relation of bulge bright- 
ness distribution i|Kormendvlll977l ). It is known that a realis- 
tic galaxy merging requires some sort of feedback (e.g. galac- 
tic winds) to prevent star formation runaway in the central 
regions. We did perform preliminary simulations including 
such winds, and find that simulated bulges obey the Kor- 



mendy relation. The problem is that these simulated bulges 
have less particles in their inner regions causing a dynamical 
relaxation time which is too short for our needs. 



7 DISCUSSION 

The standard ACDM model naturally produces massive 
classi c al bul ges in large spiral galaxies (|Scannapieco et all 
120091 . I2010I ). This is likely the result of the signifi- 
cant late violent assem bly of galactic material in ACDM 
|Peebles fc Nusserll20ldl ). Thus ACDM faces two problems 
when compared to observations of nearby galaxies. First, 
there are too many giant pure disk galaxies with none or tiny 
bulges (even pseudo-bulges). Second, there are too many 
large pseudo-bulges that are rotationally supported, in con- 
trast to the tendency of major mergers to produce classical 
bulges with little rotational support. Here we have demon- 
strated that the existence of pure disk galaxies could be 
accompanied by the existence of large pseudo-bulges as nat- 
ural outcome of mergers of pure disks. 

In contrast, the formation of pure disk dwarf galax- 
ies could be explained by Supernovae (SN) feedback. SN 
feedback ejects gas from the central region, releasing most 
central stars from the gravitation grip in the already shal- 
low gravitational potential wells of these s mall galaxies 
IIGovernato et al. 1 120101 : iNavarro et all Il996al ). The forma- 
tion of large pure disks within the ACDM model requires 
significantly more contrived conditions. Although our simu- 
lations are done with large original disk galaxies, we show 
the formation of pseudo-bulges only requires the existence 
of gas rich dwarf disk galaxies. 

Gas rich last major mergers may produce the ob- 
served abundance of pseudo-bulges in the local Universe. 
Pseudo-bulges initially created by a major merger event 
may continue to grow over long periods of time by regu- 
lar secular processes. In this picture pseudo-bulges are ex- 
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Figure 9. Percentage of solid angles giving a rotationally sup- 
ported remnant, as function of time. Black crosses, close green cir- 
cles, and closed red triangles correspond to g70qlo3n2, g70qlo3, 
and g70qlol, respectively. 



pected to contain stellar populations younger than classi- 
cal bulges - not only due to the later secular evolution, 
but also because in this scenario a significant fraction of 
classical bulges is expected to form by merging of pseudo- 
bulges. This expectation is c onsistent with observations 
(|Kormendv fc Kennicuttll2004 l. Further, according to this 
picture the fraction of classical bulges should increase with 
the density of the environment, because the galactic merg- 
ing rate depends on environment density; Higher fractions 
of mergers of pseudo-bulges into into bulges is expected in 
denser regions. Th e effect seems to be co nsistent with cur- 
rent observations (|Kormendv et al.l l2009h but further evi- 
dence is required. 

As a final note, iKormendv et all (|201ll ) confirm the cor- 
relation between black-hole (BH) mass and their host classi- 
cal bulges. However they find that BH masses correlate little 
with pseudo-bulges and not at all with their host disks. The 
common interpretation these results is that BHs co-evolve 
with their host classical bulge, but not with their host disk. 
Continuing this line of thought, the finding that BH mass 
may correlate with its host pseudo-bulge but not with its 
host disk means that some pseudo-bulges may not evolve 
from disks by secular evolution. The pseudo-bulge forma- 
tion model presented in this work may supply the mecha- 
nism needed to form such pseudo-bulges: pseudo-bulges that 
did not evolve from their host disk. 

All codes used in this paper (except for the modified 
version of Gadget2) are available at 
http : //physics . technion. ac . il/~kari/| 
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Figure Al. Equilibrium of the DM and the stellar disk (SD) of 
the gas-less galaxy gO. Mass distribution (top) and rotational ve- 
locity and velocity dispersions of the stellar disk (bottom) shown 
at t = and t = 1.05 Gyrs corresponding to llttji/n where the 
dynamical time t^ yn is computed at radius r = 6kpc. The curve 
SD ra d is the mass of the disk contained in a sphere of radius r. 
The curve SD ver t shows the stellar mass contained between two 
planes at distance r above and below the mid-plane of the disk. 



APPENDIX A: THE RECONSTRUCTION OF 
ISOLATED GALAXIES IN EQUILIBRIUM 

The simulated galaxies contain both a collision-less (DM and 
stars) and a collisional (gas) component. We initially build a 
galaxy in detailed equilibrium, in which each component is in 
kinematic equilibrium regardless of the others. We describe 
the method used for the collision-less components first, and 
then the method used for the gaseous component. 

Equilibrium initial conditions for N-body simulations 
are solutions to the collisionless Boltzmann equation. 



There are numerous methods to obtain particle posi- 
tions and velocities given some spatial or phase restric- 
tions (e.g. ISchwarzschildl jl979l: I Kuiikcn fc Dubinskil Il995l : 



I Widrow fc Dubinskill2005l ; lHernquistlll993l ). most methods 
have some severe restrictions, and are restricted to either a 
single-component, some special symmetry, collisionless com- 
ponents (do not allow for inclusion of a heavy gaseous com- 
ponent), or are merely approximate. We now describe a 
method which is efficient, general, and accurate. We de- 
velop this method mainly as a derivation of the itera- 
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Figure A2. Left panel: Equilibrium mass profiles for the DM, stellar disks (SD) and gas (G). In each panel, the curves to the left 
show the mass enclosed between two planes parallel to the stellar disk, enclosing it, and each at a distance r from it. Solid blue, dashed 
red and dotted black lines, respectively correspond to t = 0, 3td yn , and t = 6td yn - Right panel: Velocity moments equilibrium of stellar 
component, corresponding to the SD curves in the left panel. The notation of the curves is the same as in the left panel. 



tive m ethod (IM), which was introduced bv lRodionov et all 
This method belongs to the general family of 
"Made-to-Measur e " (M 2M) methods first introduced by 
ISver fc Tremaind ljl99tj ). The basic algorithm of the IM 
method is: 

(i) Generate a particle realization of all components of 
the galaxy. This realization should include only particle po- 
sitions, while velocities are not essential, although certain 
physically motivated choices would boost convergence. We 
call this realization the "original" snapshot. 

(ii) For each galactic component, use an N-body code to 
move the particles forward to a new snapshot. The integra- 
tion is done for a dynamically significant period of time, r. 
All other components are unaltered and their gravitational 
potentials are assumed static. 

(iii) Transfer velocities from the new to original snapshot. 
Assign each particle an integer number n, equals the number 
of times the particle has been used for velocity transfer. For 
every particle in the original snapshot, locate the closest n p 
particles in the new snapshot. Among those particles, find 
the nearest one with the smallest n u . Transfer the velocity of 
this particle to the original particle. This is done separately 
for every galactic component. 

(iv) repeat from step (ii) until sufficient convergence is 
achieved. 



equations (e.g. ISpringel et al.l [20051 ) . may result in compo- 
nent ^sthatare out of equilibrium in their central regions. 

iDehnenl (j2009h developed a modified M2M method 
(MM2M) in order to alleviate two main problems of the 
IM method. The first is a scale problem- for every galactic 
component, the regions which are close to the galactic centre 
have a much shorter dynamical time than the outer regions. 
When simulating with an N-body method in the IM, the in- 
ner regions achieve sufficient convergence much faster than 
the outer regions, yet they are further simulated for a long 
time until the whole model, including outer regions, achieve 
sufficient convergence. Second, The method requires a full 
N-Body simulation with time-cost comparable to the simu- 
lation itself; i.e., the method requires a computational effort 
similar to that of the si mulation. These problems are not 
straight forward to solve; iRodionov fc Orlovl (|2008l ) showed 
that using a static potential instead of integrating with the 
N-Body method results in unprovable and non-physical solu- 
tions to the Boltzmann equation, such as oscillating velocity 
profiles. 

We now describe the method used in this work, which 
we term the "Modified Iterative Method" (MIM). This 
method has equivalent benefits to MM2M, albeit it is 
achieved with a simpler technique, and is also used and 
tested with massive gaseous disks. The algorithm is as fol- 
lows: 



This method slowly assigns a developed kinematic state 
which corresponds to a given matter density. This method 
is also very flexible, and easy to implement. The kinematics 
may be constrained in many ways (e.g. isotropic or rotating 
bulges). This method has many advantages over approxi- 
mate ones to produce multi-component gravitating systems 
in equilibrium; For example, methods using a Gaussian ve- 
locity distribution with dispersion derived from the Jeans 



(i) As in the IM method, generate a particle realization 
of all components of the galaxy 

(ii) Integrate in time all particle trajectories in the global 
static potential. Integration time Ti = nitdi is individual 
per-particle, with tdi being some local dynamical time, and 
rii some small integer number common to all particles. 

(iii) Transfer velocities from new to original snapshot: this 
is done like in the IM method, except copied velocities are 
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given a "penalty" in form of addition of some random ve- 
locity proportional to the number of times it was used, n u . 



(iv) repeat from step (i) until sufficient convergence is 
achieved. 

The changes purpose are to better explore the phase-space 
to achieve better and faster conversion. With these changes, 
gravitational systems converge to a physically probable 
Boltzmann solution even when integrating orbits using a 
static potential. 

In order to include a massive gaseous disk in equilib- 
rium, we integrate the collisionless particle trajectories un- 
der the static potential of a gaseous disk. While the radial 
structure of this disk is arbitrarily chosen, the vertical struc- 
ture is a self consistent solution to the hydrodynamical and 
Poisson equations. To achieve this, we modified the N-Body 
code to allow the gaseous particle to move only in the ver- 
tical direction. In this study, gas particles are initially posi- 
tioned as a realization of eq. U Applying a dissipative term 
to their velocity, and applying the static potential from the 
collisionless components, in addition to their self-gravity, the 
particles rapidly sink to equilibrium. We find that conver- 
gence can be accelerated if using a vertical scale factor zq 
that is much smaller than the equilibrium one. This is be- 
cause the driving force in this configuration is the pressure, 
which in this case is much stronger than gravitational. 

This method allows for fast construction of galaxies 
with massive gaseous disks, is arguably simpler to imple- 
ment compared to methods with similar purpose, and more 
importantly, gives an exact solution to the Boltzmann equa- 
tion. This method is very flexible, can be used with arbitrary 
spheroidal components, and requires no analytical knowl- 
edge of their kinematic properties. 

We test the ability of the MIM at producing equilibrium 
disk galaxies. Fig. lAll shows the results of simulating the gas 
free galaxy gO in isolation for 11 dynamical times at a nom- 
inal radius of 6kpc, where td = 0.15 Gyrs. Fig. I All shows 
that there is no significant change in mass structure, or ve- 
locity structure, over this time period. In order to test the 
performance of the method when a gaseous disk is included, 
we run the additional 4 original galaxies for 6 dynamical 
times, as measured at our nominal radius of 6kpc. The re- 
sults are shown in Fig. IA21 where the solid blue, dotted 
black, and dashed red curves correspond to t — 0, t — Stdyn, 
and t = 6td yn - When comparing to the gO simulation, it is 
evident that the massive gaseous disk damages the quality 
of equilibrium achieved with the MIM method. The reason 
behind this is that gas instabilities change the structure of 
the gaseous disk (as seen by its the changing thickness). 
However, equilibrium is reached within a short time, as seen 
by the small difference between the dashed red lines and the 
dotted black lines. This result applies both for mass and ve- 
locity structure of the galaxies, as shown in the right panel 
of the same figure. In the merging simulations the galaxies 
are initially separated far enough as to allow the gaseous 
disks to achieve better equilibrium. 



